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SUMMARY 


The effect of interaction between combustion processes and structural 
deformation of solid propellant was considered. The combustion analysis was 
performed on the basis of deformed crack geometry, which was determined from 
the structural analysis. On the other hand, input data for the structural 
analysis, such as pressure distribution along the crack boundary and ablation 
velocity of the crack, were determined from the combustion analysis. The inter- 
action analysis was conducted by combining two computer codes, a combustion 
analysis code and a general purpose finite element structural analysis code. 


INTRODUCTION 


In recent years, much attention has been focused on the investigation of 
the coupling effect between combustion phenomenon and mechanical behavior of solid 
propellant. The solution of problems of this type can further better understand- 
ing of the transient combustion processes inside solid propellant cracks, which 
may significantly affect the performance of a rocket motor. The combustion pheno- 
menon inside the crack of solid propellant is strongly influenced by the crack 
geometry as the material is being deformed and burned away. Generally, there are 
two major reasons for alteration of the crack geometry: 1) mass loss due to 

gasification of propellant surface along the crack during the combustion process, 
and 2) mechanical deformation of the propellant due to pressure. 

On one hand, both the burning rate and mechanical deformation are governed 
by pressure acting on the crack surface. On the other hand, a change in crack 
size will cause the pressure distribution to vary. The pressure distribution 
will strongly influence the deformation and stress concentration at the crack 
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tip, which in turn will affect the manner of the crack propagation. It is, 
therefore, apparent that the pressure distribution and the change in crack geo- 
metry are strongly interdependent. 

In the past, combustion and structural analyses of solid propellant were 
conducted independently, with the result that interaction effects were completely 
ignored. As noted above, such interaction effects can be quite important, especi- 
ally when the deformation is large as compared to the original crack-gap width. 

The deformation response of the material is categorized as linearly viscoelastic. 

It is, therefore, the intent of this paper to present a method of analysis for the 
combust ion- structural interaction in a linear viscoelastic medium. To this end, 
three major tasks are involved: 1) combustion analysis to model the transient 

combustion process, 2) viscoelastic analysis in conjunction with moving boundary, 
and 3) linkage of the two analyses. 

For the combustion analysis, investigations of certain aspects of combustion 
processes have been made. Taylor [1] conducted experimental tests to study the 
convective burning of porous propellants with closed- and open-end boundary condi- 
tions. Belyaev et al. [2] showed that the burning of propellant inside a narrow 
pore may lead to an excess pressure buildup. In a later study, Belyaev et al. [3] 
made a series of experimental tests to determine the dependence of flame-spreading 
rate on crack geometry, propellant properties, boundary conditions, and combustion 
chamber pressures. Cherepanov [4] stated that as a result of the impeded gas flow 
in a sufficiently narrow and long cavity, the pressure reaches such high values 
that the system becomes unstable. From his work, Godai [5] indicated that there is 
a threshold diameter or critical width of a uniform cavity below which flame will 
not propagate into the crack. Krasnov et al. [6] investigated the rate of pene- 
tration of combustion into the pores of an explosive charge. Jacobs et al. [7,8] 
studied the pressure distribution in burning cracks that simulate the debonding 
of solid propellant from the motor casing. 

Although results of previous experiments were of interest, no sound theore- 
tical model was developed. In this study, a theoretical model was established 
for predicting rate of flame propagation, pressure distribution, and pressuriza- 
tion rate inside the crack. Two sets of coupled partial differential equations 
were obtained: one from mass, momentum, and energy conservation of the gas phase 

of the propellant product in the void region adjacent to the crack surface; the 
other from consideration of solid-phase heat conduction. Due to the mathematical 
complexity of governing equations and boundary conditions involved, the finite 
difference method was used to obtain the solution for the combustion analysis. 

In the numerical solution, the boundary conditions, which vary with time, are 
specified in terms of the changing crack geometry, which in turn is found from the 
structural analysis. In addition, the pressure distribution along the crack sur- 
face, varying as a function of time, was obtained from the analysis and was used 
as input for the structural analysis. 

For structural analysis, different approaches have been taken previously 
in solving (analytically or numerically) several moving boundary problems in linear 
viscoelasticity. Lee et al. [9] obtained a solution for the pressurization of 
an annihilating viscoelastic cylinder contained by an elastic casing in which 
the material was assumed to be a Kelvin model in shear and incompressible in bulk. 
Arenz et al. [10] performed a similar analysis for a sphere. Corneliussen et 
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al. [11,12] presented solutions for a spinning, annihilating, viscoelastic 
cylinder with free outer boundary. Since the constraining case was not included 
in their analyses, the stress distribution was independent of the material pro- 
perties. With the assumption of a standard linear solid model, Shinozuka [13] 
presented the analytical solution for a case-bonded pressurized viscoelastic 
cylinder. More generalized solutions were obtained by Rogers et al. [14] for a 
class of linear viscoelastic problems by using the numerical integration scheme. 
Schapery [15] also developed a general method for solving moving boundary problems. 
In his approach, the moving boundary condition was replaced by a fictitious non- 
moving boundary subjected to a time-dependent pressure. Later, Christensen et 
al. [16] obtained a series solution for the stresses of the same problem. As 
noted above, most of the analytical solutions were available for viscoelastic prob- 
lems of simple geometry. For complex geometry, the finite element method has 
proven to be most useful - 

Application of the finite element method for solving viscoelastic problems 
is not new; reports of such work can be found, for instance, in references [17-20]. 
However, most of the previous work did not consider the effect of moving boundary, 
an important feature for the structural analysis of solid propellant. Sankaran 
and Jana [21 ] presented a technique for the solving of axisymmetric viscoelastic 
solids with moving boundary. In their approach, the finite element mesh corres- 
ponding to the new boundary was re-generated, while the stress-strain histories 
and material properties were assumed to be carried over from those of the previous 
time increment. This assumption is valid only if the time increment is very small. 
An algorithm for automatically tracking ablating boundaries was given by Weeks and 
Cost [22]. All previous work dealing with moving boundary viscoelastic problems 
lacks both the appropriate treatment of material properties, and stress-strain 
histories for the newly generated mesh. It is the purpose of this paper to pre- 
sent such a treatment. 

Three major features must be included in the structural analysis for a solid 
propellant: 1) proper modeling of viscoelastic behavior, 2) tracking of 

ablating boundary in order to generate new finite element meshes, and 3) treat- 
ment of the material responses (i.e., stress-strain histories and material pro- 
perties) for the new mesh. All of these features have been incorporated into a 
nonlinear finite element program called NFAP [23]. Combustion and structural 
programs were combined in order to make possible an interaction analysis. Numeri- 
cal results are presented to demonstrate the effect of interaction between 
combustion and structural responses of the material. 


COMBUSTION ANALYSIS 


The theoretical model was developed to simulate the combustion phenomenon 
inside a propellant crack, which is located in a transverse direction to the 
main flow of the rock chamber. During the course of derivation, the following 
assumptions are made: 

1) All chemical reactions occur near the propellant crack surface, and the 
combustion zone is so thin that it is considered a plane. 


2) Rate processes at the propellant surface are quasi-steady in the sense that 
characteristic times associated with the gaseous flame and preheated pro- 
pellant are short in comparison to that of pressure transient variation. 
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3) Gases flowing in the propellant crack obey the Clausius or Noble-Abel 
equation of state. 

4) Bulk flow in the pore is one-dimensional [24]. 

To describe gas-phase behavior inside a solid propellant crack, mass, 
momentum, and energy equations in unsteady, quasi-one -dimensional forms have 
been developed, based upon the balance of fluxes in a control volume within the 
propellant crack. 

The mass conservation equation is 
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The energy conservation equation written in terms of the total stored energy 
(internal and kinetic) per unit mass, is 
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The conservation equations are further simplified by an order of magnitude 
analysis in which the following terms are negligible: 1) forces between mole- 

cules due to viscous normal stress in axial direction; 2) viscous dissipation 
and rate of work done by the force caused by viscous normal stresses in the energy 
equation; and 3) axial heat conduction between gas molecules in the energy 
equat ion • 

The propellant surface temperature at a fixed location along the crack before 
the attainment of ignition is calculated from the solid-phase heat conduction 
equation written in unsteady one-dimensional form: 
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where the length variable y is measured perpendicular to the local propellant 
crack surface. Initial and boundary conditions are 
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The heat conduction equation is solved by using an integral method [25] 
which employs a third-order polynomial, or by direct numerical solution of Eqs. 
(4-7) with variable mesh size in the subsurface. 


For the gas phase, the Noble-Abel equation is used for the equation of 
state: 


P^^ - b) = RT (8) 

The gas-phase equations, i.e. Eqs. (1), (2) and (3), are non-linear, inhomo- 
geneous, partial differential equations. Along with the partial differential 
equation for the solid phase (Eq. (4)), they are solved simultaneously, using the 
finite difference method. The derivation described above was implemented into a 
computer program, crack combustion code (CCC) by Kuo et al- [26]. 


STRUCTURAL ANALYSIS 


To conduct the structural analysis of the solid propellant, three main 
features must be included in the numerical formulations: 1) modeling of 

viscoelastic material behavior, 2) simulation of ablating boundary, and 
3) treatment of material responses by an interpolation scheme. Each feature 
is outlined below. 


Viscoelastic Material Model 

The material behavior of the solid propellant is assumed to be visco- 
elastic in shear and elastic in bulk. Only the isothermal condition is considered. 
The stress-strain relations with zero initial conditions are written in two parts. 

1) Shear behavior: 
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where Is the relaxation modulus In shear. For most viscoelastic 
materials, it is usually considered 
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2) Bulk behavior: 
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As discussed in [27], the incremental stress-strain relations in matrix form 
are written as 
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and the term C.. has a recursive relationship, i.e., 
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The advantage of Eq. (19) is that all of the strain history can be obtained by 
referring only to information in the previous time step, thus reducing computer 
storage and numerical calculations. 


From the virtual work principle and the relationship of Eq. (12), the finite 
element equilibrium equations for a typical time interval [t, t + At] can be 
derived as 
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Simulation of Ablating Boundary and Mesh Generation 

Burning of the propellant causes a significant change in geometry, thus 
presenting complications in finite element structural analysis. The effect of 
ablating boundary is accounted for by redefining the finite element mesh at speci- 
fied time intervals. This involves two stages of calculations: 1) tracking 

of the ablating boundary, and 2) generation of new finite element mesh. With 
some modifications, the procedures adopted herein are similar to those presented 
in [22]. 

Consider a structural geometry with ablating boundary. The spatial posi- 
tions of the ablating boundary are determined by the ablation velocities which 
are found from the combustion analysis at discrete times. It is assumed that 
the ablation occurs always in the direct on normal to the boundary. For struc- 
tural analysis, the entire surface is divided into an ablating part and a non- 
ablating part; each part is formed by discrete line segments joining at the nodes 
of the finite element mesh. The new position of each line segment is located from 
the given ablating velocity. Consequently, the new boundary nodes are determined 
by calculating the intersections of two subsequent new line segments. Likewise, 
the nodes at the intersections of new ablating and non-ablating boundaries are 
then determined. 

During the locating process, however, some of the boundary nodes may not lie 
on the new boundary and thus must be eliminated. If the distance from the tip 
of the normal vector at a new nodal position to any node on the original boundary 
is less than the value of the normal itself, the node is removed. 
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In general, the total number of nodes on the boundary at discrete times 
will be different because some of the nodeshave been removed. However, in 
the analysis it is more convenient to generate a finite element mesh similar to 
the original one so that interpolation of material response can be made. One 
way to accomplish this is by keeping the number of boundary nodes constant. 
Consequently, the boundary nodes are redistributed between two discontinuity 
points which are specified in the input data in such a way that the lengths 
of the new line segments have the same ratio as those of the original lines. 

Once the new boundary nodes are defined, an automatic mesh generation 
scheme is used to create the interior nodes for further analysis. Because of 
its flexibility in obtaining a desirable mesh, a Lap lacian- isoparametric grid 
generation scheme [28] is utilized. However, this method is limited to a 
geometry bounded by four sides. A finite element mesh is shown in Fig. 1. The 
coordinates of the i-th interior node can be expressed in terms of those of 
neighboring nodes by 
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of interior nodes. 


Setting up the equations for each interior node yields two systems of simul- 
taneous equations. It is observed that the resulting systems of equations are 
banded and symmetric. The Gaussian elimination scheme is employed to solve for 
the coordinates of the interior nodes. 


Interpolation of Material Responses 

As seen from Eq. (16), the stress increment Aa for the time interval 
[t, t+At] varies with material properties and with the strain history at both 
current and previous time steps. When the region of an element changes over a 
period of time due to ablation, the material response history of the new elements 
is lost and must be determined by an interpolation procedure from the old ele- 
ments at previous time steps. Accordingly, the interpolation procedure is carried 
out on the element level. For calculations, the material responses are separated 
into two groups : the first includes such variables evaluated at the Gaussian 

integration points, i.e., Aa.., Ae.. and C?.; the second includes the nodal 
displacements which are evaluated ai nodaT points. In the present calculations, 
two limitations are imposed: 1) eight -node quadrilateral elements are used 

throughout the analysis; and 2) the four sides of each element remain straight 
before and after ablation. 


1) Interpolation of Gaussian variables - It is noted that the quadratic dis- 
placement approximation of an eight -node element yields a linear strain varia- 
tion. With this fact in mind, the quantities of Gaussian variables at nodal 
points are first evaluated for every old element. As shown in Fig. 2a, b, this 
can be done by using the linear isoparametric shape functions, namely. 
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For each new element, the local coordinate^ (r^, s) of the k-th Gaussian 
point are known. The global coordinates, (y, , z ), of that point are, therefore, 
computed by using the following equations: 
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After (y, j z, ) are found, the old element to which the point belongs must be 
identified. search process based upon the values of r’ and s’ is developed for 

this purpose. The search starts from the old element which corresponds to the 
neighboring elements . Equations for such calculations are given by 
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Fig. 2c shows how to identify the element to which the points, (r’, s’), belong. 
Once the location of the point is verified, an interpolation procedure is per- 
formed, using the relationship 

_ 4 

w, = ^ h.(r’, s’) * w.' (29) 
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2) Interpolation of nodal displacements - A similar procedure to that explained 
above is also used to determine the position of the node in question with refer- 
ence to the old element. However, the interpolation procedure in Eq. (26) is no 
longer necessary since the nodal displacements are known. The nodal displacements 
of the new mesh are computed from 
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where h. are the standard quadratic isoparametric shape functions. 


75 


All formxilations discussed in this section have been implemented into a 
general purpose nonlinear finite element program called NFAP for conducting 
viscoelastic analysis of solid propellant with ablating boundary. Some numer- 
ical examples are presented in a later section. 


COUPLING EFFECT 


For structural analysis, the boundary condition along the crack geometry 
is defined by pressure distribution which varies with time, and ablation velocity; 
both are determined from combustion analysis. In the combustion analysis, the 
regression rate of the propellant is dependent on the deformed crack geometry. 
Therefore, the two processes are strongly interdependent. Such a coupling effect 
is obtained by combining the analysis of two computer programs: a crack combus- 

tion code (CCC) and a structural analysis code (NFAP) . Both codes were developed 
independently to facilitate program verifications. Linkage of the two codes was 
made subsequently. 

The coupling effect considered in the present analysis is limited to the 
major parameters, namely pressure loading, ablation velocity, and crack deforma- 
tion. Pressure and ablation velocity are calculated by the CCC at each nodal point 
located on a one-dimensional grid along the length of the crack. The analysis 
of crack combustion incorporates the crack geometry variation caused by both 
mechanical deformation and mass loss through gasification of the propellant 
surface. Once the gas-phase equations are solved and the pressures and ablation 
velocities along the crack are calculated for a given time t, the data are trans- 
ferred to the NFAP as the input information. NFAP then simulates the updated crack 
geometry from the ablation velocities and generates a new finite element mesh. 

With the new mesh and pressure data, NFAP updates the stiffness matrix and inter- 
polates material responses for conducting a quasi-static analysis at time t. After 
obtaining the deformation, the change in the crack width at each finite different 
node is calculated and added to the existing crack width. Since the crack width is 
the input of the combustion analysis, one cycle of calculations is thus completed. 
The same procedure is followed for every specified time increment. 


EXAMPLES 


For program verification and demonstration of its analysis capability, three 
sample problems were run either by NFAP alone or in the combined NFAP/ CCC program. 
The results of the analysis are discussed in the following. The nxjmerical results 
obtained from CCC alone are contained in reference [26]. 

1. A Reinforced Thick-walled Cylinder 

Figure 3 shows a cylinder of viscoelastic material bonded by a steel casing 
and subjected to a step- function internal pressure. The example was selected 
because it is composed of two different materials and the analytical results 
are readily available for comparison. Only five eight-node axisymmetric elements 
were used to model the cylinder. The material properties of the elastic casing 
are 

E = 2.068 X 10^ MPa V = 0.3015 
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The material properties of the viscoelastic core are defined by 


K = 689-5 MPa = 51.71 * exp (-0.lt) MPa 

In Fig. 4, the variations of circumferential stresses with time are plotted 
for comparison with the analytical solution obtained in [9]- It is observed that 
both solutions agree very closely. This problem was analyzed previously by 
Zienkiewicz et al. [18], using strain rate formulation of the finite element 
method. However, the formulation presented in the present paper is more easily 
Incorporated into the NFAP program. 

2. A Star-shaped Solid Rocket Motor 

As an application of the present approach in dealing with the moving boundary, 
a star-shaped solid rocket motor was analyzed by assuming both a constant and 
ablating inner boundary. The configuration and finite element mesh are shown 
in Fig. 5, and the material properties of outer casing and inner propellant are 
identical to those of the first example. Taking advantage of the symmetry 
condition, only a 30°-sector was modeled by finite element mesh. The contours 
of maximum compressive stress analyzed by constant inner boundary at various 
times are shown in Fig. 5. Comparing the present results with those of [18], it 
is evident that the general pattern is quite similar but that some small differences 
do exist. Since the geometry of the rocket motor in [18] was not clearly defined, 
the difference in dimension used in these two analyses could be the cause of such 
deviations. 

The actual case of a solid rocket motor can be modeled more closely by con- 
sidering the inner boundary being ablated. Figure 6 shows the contours of maximum 
compressive stress predicted by NFAP, using the option of moving boundary. The 
results obtained are quite different from those of [18]. However, observing the 
differences between Figs. 5 and 6, we can conclude that the results obtained by 
NFAP are quite reasonable. The solution reveals that the high stress region 
obtained for ablating boundary propagates faster than that with non-ablating 
boundary . 

3. A Propellant Crack Specimen 

As a final example, a propellant crack sample was analyzed, using the 
combined NFAP/CCC program to demonstrate the coupling effect. The initial geo- 
metry and finite element mesh generated by NFAP is given in Fig. 7. The crack 
is 0.15 m long and the initial gap-width is 0.89 mm. The web thickness is 
8 mm along the crack and 20 mm at the tip- Because of symmetry, only half of the 
sample was modeled by 80 plane strain elements. The shear relaxation modulus of 
the propellant was assumed to be 

G^(t) = 1.461 + 7.43 * exp (-.095t) MPa; and K = 4,826 MPa. 

Calculated pressure distributions at various times, from the CCC alone, are 
given in Fig. 7. The burning phenomenon of the propellant can be briefly des- 
cribed as follows. The pressure in the chamber increases with time, causing 
the hot gases to penetrate further into the crack. As time passes, the pressure 
wave travels along the crack and is reflected from the closed end. At about 
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200 ys, the pressure front has already reached the tip and is reflected, causing 
pressure at the tip to increase. 

Figure 8 shows the results obtained from the combined NFAP/CCC program- 
During the initial period, the general trend of the pressure distribution is simi- 
lar to that from convective burning analysis alone. However, as time progresses, 
noticeable differences between the two cases begin to appear. Up to 200 ys, the 
pressures obtained from the combined analysis are lower, except near the crack 
entrance region. At t = 300 ys, two pressure peaks appear. At t =» 325 ys, three 
pressure peaks appear. These pressure peaks are caused by the partial closure of 
the gap. The deformation pattern of the propellant is quite irregular because 
of the uneven distribution of the pressure along the crack surface. The elements 
at the crack entrance are compressed by the high chamber pressure, which results 
in the propellant being pushed into the crack. Since chamber pressure increases 
more quickly than pressure inside the crack, the propellant is pushed toward the 
lower pressure region inside the crack. The mechanical deformation of the 
propellant causes narrowing of the crack width, and consequently results in a local 
crack closure. This local gap closure manifests itself in a pressure peak. The 
localized pressure peaks or gap closures move along the crack. At t = 325 ys, this 
localized pressure phenomenon becomes evident at x/L = 0.167, 0.433, and 0.633. 


CONCLUSION 


The computer program for evaluating the coupling effect between convective 
burning and structural deformation was developed by combining the Crack Combustion 
Code and a Nonlinear Finite-Element Analysis Program. In structural analysis, 
the linear viscoelastic material model, together with the capabilities of simu- 
lating ablating boundary and interpolating material responses, was considered. 
Also, the coupling effect estimated by the combined analysis shows some signifi- 
cant interaction between the combustion and mechanical deformation. This pheno- 
menon will be verified further by future experiments. 


78 



SYMBOLS 


1. Combustion Analysis 

A = cross-sectional area of crack 
P 

= body force 
b = co-volume 

c^ = specific heat at constant pressure 
E = total stored energy 

h = local convective heat-transf er coefficient 
c 

local convective heat-transfer coefficient over propellant surface 

h^^= local convective heat-transfer coefficient over nonpropellant port wall 

h^ = enthalpy of combustion gas at adiabatic flame temperature 

p, = burning perimeter 
ri 

P = wetted perimeter of port 
w 

p = static pressure 

R = specific gas constant for combustion gases 

r, = burning rate of solid propellant, including erosive burning contribution 

D 

X = temperature (without subscript, static gas temperature) 

= adiabatic flame temperature of solid propellant 

X .= initial propellant temperature 
pi 

X = propellant surface temperature 
ps 

X^^= nonpropellant wall surface temperature 
t = time 
u = gas velocity 

velocity of propellant gas at burning surface 
X = axial distance from propellant crack opening 
y = perpendicular distance from propellant surface into solid 
a = thermal diffusivity 
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Y = ratio of specific heats 
X = thermal conductivity 
U = gas viscosity 

p = density (without subscript, gas density) 

T = shear stress on port wall 
w 

T = normal viscous stress 

XX 

0 = angle measure, in a counterclockwise direction, at lower side of 

propellant , degree 

Subscripts 

1 = initial value 

pr = solid propellant (condensed phase) 
c = rocket chamber 
2. Structural Analysis 

stress deviators 
strain deviators 
stress tensor 
shear relaxation modulus 
bulk modulus 
a quantity at time t 
incremental stress 
incremental strain 

equivalent initial stress vector due to viscoelastic behavior 
stiffness matrix 
transpose of matrix 

number of terms of series in relaxation modulus 


S. . « 

= 

K 

( )^ = 

{Aa} = 

{Ae} = 

{a } = 
o 

[K] = 

C = 

M 
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material constants in relaxation modulus 

{Av} = increment of nodal displacement vector 

[Dve^ = viscoelastic material matrix 

Cb] — strain-nodal displacement transformation matrix 

= values of Gaussian variables at k-th integration point referred 
to old element 

= values of Gaussian variables at 1-th nodal point referred to 
old element 

(r ’,s ’) = local coordinates of k-th integration point referred to old 
element 

(y^,z^) = global coordinates of 1-th nodal point referred to new element 
= i-th nodal displacement referred to new element 
= i-th nodal displacement referred to old element 
(r*,s’) = local coordinates of point in equation referred to old element 


Sin 


81 



REFERENCES 


1. Taylor, J,W. : The Burning of Secondary Explosive Powders by a Convective 

Mechanism. Trans. Farad. Soc. 58, 1962, p. 561. 

2. Belyaev, A.F. ; Korotkov, A. I.; Sulimov, A.A. ; Sukoyan, M.K. ; and Obemin, 

A.V.: Development of Coribustion in an Isolation Pore. Combustion, Explosion 

and Shock Waves, Vol. 5, Jan. -March 1969, pp. 4-9. 

3. Belyaev, A.F.; Bobolev, V.K. ; Korotkov, A.A. ; Sulimov, A.A.; and Chuiko, 

S.V.: Development of Burning in a Single Pore, Transition of Combustion of 

Condensed Systems to Detonation, Chap. 5, Pt.A. ,. Science Publisher, 1973, 
pp. 115-134. 

.4. Cherepanov, G.P. : Combustion in Narrow Cavities. J. of Appl. Mech., Vol. 11, 

1970, pp. 276-281. 

5. Godai, T. : Flame Propagation into the Crack of a Solid-Propellant Grain. 

AIAA Journal, Vol. 8, July 1970, pp. 1322-1327. 

6. Krasnov, Yu. K. : Margulis, V.M. ; Margolin, A.D.; and Pokhil, P.F. : Rate 

of Penetration of Combustion into the Pores of an Explosive Charge. Combus- 
tion, Explosion, and Shock Waves, Vol- 6, July-Sept. 1970, pp. 262-265. 

7. Jacobs, H.R. ; Williams, M.L.; and Tuft, D.B.: An Experimental Study of the 

Pressure Distribution in Burning Flaws in Solid Propellant Grains. Univ. of 
Utah, Salt Lake City, Utah, Final Report to Air Force Rocket Propulsion 
Laboratory, AFRPL-TR- 72-108 , UTEC DO 72-130, Oct. 1972. 

8. Jacobs. H.R. ; Hufferd, W.L. ; and Williams, M.L. : Further Studies of the 

Critical Nature of Cracks in Solid Propellant Grains, AFRPL-TR- 7 5-1 4, 

March 1975. 

9. Lee, E.H.; Radok, J.R.M. ; and Woodward, W.B.: Stress Analysis for Linear 

Viscoelastic Material. Trans. Soc. Rheol., Vol. 3, 1959, pp. 41-59. 

10. Arenz, R.J.; and Williams, M.L. : The Stresses in an Elastically Reinforced 

Pressurized Viscoelastic Sphere with an Eroding Boundary. Proceedings 
20th Meeting Joint Army-Navy-Air Force Physical Properties Panel, Johns 
Hopkins Univ. , Baltimore, Md., 1961, p. 143. 

11. Comeliussen, A.H.; and Lee, E.H.: Stress Distribution Analysis for Linear 

Viscoelastic. Materials. International Union of Theoretical and Applied 
Mechanics Symposium on Creep, 1960, pp . 1-20. 

12. Comeliussen, A.H, ; Kamowitz, E.F. ; Lee, E.H. ; and Radok, J.R.M. : Visco- 

elastic Stress Analysis of a Spinning Hollow Circular Cylinder with an 
Ablating Pressurized Cavity. Trans Soc. Rheol., Vol. 7, 1963, pp. 357-390. 

13. Shinozuka, M. : Stresses in a Linear Incompressible Viscoelastic Cylinder with 

Moving Inner Boundary. J. Appl. Mech., Vol- 13, 1963, pp. 335-341. 

14. Rogers, T.G.; and Lee, E.H. : The Cylinder Problem in Viscoelastic Stress 

Analysis. Quart. Appl. Math., Vol. 22, 1964, pp. 117-131. 


82 


15. Schapery, R.A. : An Approximate Method of Stress Analysis for a Large Class 

of Problems in Viscoelasticity. Purdue Univ. Rept. A. and ES62-18, 1963. 

16. Christensen, R.M. ; and Schreiner, R.N. : Response to Pressurization of a 

Viscoelastic Cylinder with an Eroding Internal Boundary, AIAA Journal, 

Vol. 3, 1965, pp. 1451-1455. 

17. Chang,. T.Y. : Approximate Solutions in Linear Viscoelasticity. Ph.D. 

Dissertation, Dept, of Civil Eng., Univ. of California, Berkeley, 1966. 

18. Zienkiewicz, O.C.; Watson, M. ; and King, I.P.: A Numerical Method of 

Viscoelastic Stress Analysis. Int. J. Mech. Sci., Vol. 10, 1968, pp. 807-827. 

19. White, J.L.: Finite Element in Linear Viscoelasticity. Proceedings of 

2nd Conference on Matrix Method in Structural Mechanics, Airforce Flight 
Dynamics Lab, Wright Patterson AFB, Ohio, AFFDL-TR-68-150 , 1968, pp. 489-516. 

20. Gupta, K.K. ; and Heer, E. : Viscoelastic Structures. Structural Mechanics 

Computer Programs, Univ. Press of Virginia, Charlottesville, 1974, pp. 207-225. 

21. Sankaran, G.V. ; and Jana, M.K. : Thermoviscoelastic Analysis of Axisymmetric 

Solid Propellant Grains. J. Spacecraft, Vol. 13, 1976, pp. 641-642. 

22. Weeks, G.E.; and Cost, T.L.: An Algorithm for Automatically Tracking Ablat- 

ing Boundary. Int. J. Num. Method. Engr., Vol. 14, 1979, pp. 441-449. 

23. Chang, T.Y. ; and Prachuktam, S.: NFAP - A Nonlinear Finite-Element Analysis 

Program. Dept, of Civil Engr., Univ. of Akron, Rept. No. SE76-3, 1976. 

24. Kuo, K.K.; Chen, A.T.; and Davis, T.R. : Transient Flame Spreading and 

Combustion Process Inside a Solid Propellant Crack. AIAA Paper 77-14, AIAA 
15th Aerospace Science Meeting, Jan. 1977. 

25. Goodman, T.R.: Application of Integral Methods to Transient Nonlinear 
Heat Transfer. Advances in Heat Transfer, Vol. 1, Academic Press, New York, 
1964, pp. 51-122. 

26. Kuo, K.K. ; Chen, A.T. ; and Davis, T.R.: Convective Burning in Solid-Propellant 

Cracks. AIAA Journal, Vol. 16, June 1978, pp . 600-607. 

27. Chang, J.P.: Finite Element Analysis of Linear Viscoelastic Solids. M.S. 

Thesis, Dept, of Civil Engr., Univ. of Akron, 1980. 

28. Herrmann, L.R.: Lap lacian— Isoparametric Grid Generation Scheme. J. of 

Engr. Mech. Div., Proc. of the ASCE, Vol. 102, Oct. 1976, pp. 749-756. 


83 



Figure 2.- Interpolation scheme. 
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Figure 3.- Finite element mesh of a reinforced thick-walled cylinder 











Figure 6.- Contours of maximum principal compressive stress of 
a star-shaped rocket motor with moving boundary. 
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Figure 8.- Calculated pressure distributions for various 
times from the combined crack combustion and 
non-linear finite- element analysis program. 
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